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Abstract. A program is outlined, and first results described, in which fully three-dimensional, 
time dependent simulations of hydrodynamic turbulence are used as a basis for theoretical inves- 
tigation of the physics of turbulence. The inadequacy of the treatment of turbulent convection 
as a diffusive process is indicated. A generalization to rotation and magnetohydrodynamics is 
indicated, as are connections to simulations of 3D stellar atmospheres. 
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1. Introduction 



' John von Neumann (von Neumann 1948) proposed a way to deal with the intractable 
O ^, problem of hydrodynamic turbulence, by (1) using numerical simulation on computers 
1^ to construct turbulent solutions of the hydrodynamic equations, (2) building intuition 
from study of these solutions, and (3) constructing analytic theory to describe them. 
^ He proposed that iterating this procedure could lead to a practical understanding of 
, turbulent flow. The computer power available at that time was totally inadequate to 
compute hydrodynamics on sufficiently refined grids to produce turbulent flow; numerical 
I , viscosity restricts the effective Reynolds number. Today, computing power is adequate for 
J> ' the simulations of truly turbulent, three-dimensional (3D), time dependent, compressible 
OO , flows, so we have begun a program based upon von Neuman's proposal. 

Turbulent flow in its many guises (e.g., convection, overshooting, shear mixing, semi- 
[ convection, etc.) is probably the weakest aspect of our theoretical description of stars 
• ' (and accretion disks). The full problem that faces us includes rotation, magnetic fields, 
, and multi-fluids (to account for compositional heterogeneity, diffusion, radiative levita- 
' tion, and nuclear burning). In this paper we describe the progress made toward von 
' Neumann's goal. We plan to replace the venerable mixing-length theory (MLT) with a 
^ physics-based mathematical theory which can be tested by refined simulations and ter- 
. ^ , restrial experiment (e.g., laboratory fiuid experiments, meteorological and oceanographic 
' observations). Particularly relevant are high-energy density plasma (HEDP) experiments, 
^ . which now can access regions of temperature and density that overlap stellar conditions 
■ - - ' up to helium burning (Remington et al. 2000, 2006, Drake 2006), and deal with plasma 
and magnetic fields, just like star matter, not with an unionized fluid like air or water. 



2. Inadequacy of the Diffusion Model of Convection 

It is numerically convenient to replace convective mixing in a stellar evolution code 
by a diffusion algorithm, but this is not physically correct. The correct equation for the 
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change of composition Yi is (Arnett 1996), 

ay, + V- vy^ = -u- vr, + 7^„ (2.1) 

where the term on the left-hand side is the Lagrangian time derivative of the composition 
in a comoving spherical shell with velocity v, the first term on the right-hand side is the 

mixing due to rotation and turbulent velocities u across the Lagrangian shell bound- 
ary, and the last term is the composition change due to nuclear reactions which change 
species i. Thus, 

Tli = -YiAi + YjAj 
-YiYkpMA{(Tv) + YiY^pMa{(tv) + ■■■ (2.2) 

where the terms on the right-hand side represent all the ways in which species i can be 
made or destroyed. The advection operator 

-u • VFi (2.3) 

involves a velocity field u which is determined non-locally and a first order spatial gradient 
Vli. This is replaced by 

^|(4„V)"i.(||)l, (2.4) 

which has a second order derivative in space and a phenomenological local diffusion 
coefficient D. Except for contrived cases, these are the same (zero) only in the limit that 
composition is homogeneous. We need a major community effort to base stellar mixing 
algorithms on physics, comparable to the efforts led by Willy Fowler for nuclear reaction 
rates, so that both the advection and reaction terms are reliable. 



3. The Simulation Step 

We have simulated turbulent flow resulting from shell oxygen burning in a presuper- 
nova star. Because of the fast thermalization time (unlike the solar convection zone, for 
example) we can simulate the entire convective depth as well as the stable boundaries. 
This is a "convection in a box" approach, implicit large eddy simulation (ILES). Using 
a monotonicty preserving treatment of shocks (like PPM, see Boris 2007 and Woodward 
2007) insures that the turbulent energy moves from large scales to small in a way close to 
that envisaged by Kolmogorov 1941, 1962. Because the rate of the cascade of turbulent 
flow from large scales to small is set by the largest scales, there is no need to resolve the 
smallest scales, which are far below our grid resolution. This would not be the case if the 
nuclear burning time were shorter than the turnover time, instead of a thousand times 
longer. A detonation or deflagration, in which the turnover time is much longer than the 
reaction time, is a more difflcult problem. 

The aspect ratio is chosen to be large enough so that it has little effect on the simu- 
lation. The initial state is mapped from a ID model with sufScient care so that there is 
very little transient "jitter". The convection develops from numerical roundoff noise or 
from low amplitude seed noise. A quasi-steady state, in an average sense, develops in one 
turnover time, so that memory of initial errors is quickly lost. 

The simulations show that this oxygen shell burning is unstable to nuclearly-encrgized 
pulsations (primarily radial), which couple to the turbulent convective flow. The con- 
vective kinetic energy shows a series of pulses with an amplitude change of order of a 
factor of two. These disappear if the burning is artificially turned off. For more detail, 
see Meakin & Arnett (2007b). Meakin & Arnett (2006) find that 2D simulations which 
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include multiple burning shells show interactions between the shells; 3D simulations of 
multiple shells are planned. Neither have the pulsations, nor the interaction of burning 
shells, been included in any ID progenitor models to date. 

Another novel feature found by Meakin & Arnett (2007b) is entrainment at convection 
boundaries. The physics of the process is interesting; it involves the erosion of a stably- 
stratified layer by a turbulent velocity field, mediated by nonlinear g-mode waves (Meakin 
& Arnett 2007a). 

While these simulations do not contain an entire star, and thus limit the accuracy of the 
description of low-order modes, whole star simulations are developing enough resolution 
to exhibit turbulent flows (Brun 2009). Since we find that even modest resolution will 
give reliable average quantities (see below), we expect the "simulation step" to be soon 
generalized and extended to include rotation and magnetic fields. 

4. The Analysis Step 

The pressure, density and velocity were subjected to a Reynolds decomposition, in 
which average properties and fluctuating properties are separated. For example, for pres- 
sure, P = Pq + P' , so that averages give (P') = and (P) = Pq. Note that in general 
((P')^) 7^ 0. We use two levels of averaging: one over solid angle (the extent of our grid 
in and 0), and one over time (two turnover times). The resulting averaged properties 
have a robust behavior that was insensitive to grid size, aspect ratio, and limits to the 
size of the averaging dimension (provided it was large enough; two turnover times and 
60 degrees worked fine). 

Arnett, Meakin, & Young (2009a) find that the velocity scale is well estimated by 
equating the increase in kinetic energy due to buoyant acceleration to the decrease due 
to turbulent damping in the Kolmogorov cascade. This implies that it will be possible to 
make quantitatively correct estimates of wave generation and entrainment at convective 
boundaries. Arnett, Meakin, & Young (2009a) have shown that, in the solar case, the 
velocity scale is significantly larger than estimated by MLT (a factor of ~ 2), but agrees 
with both 3D atmospheres (Asplund 2005, Nordlund & Stein 2000 and Stein & Nordlund 
1998) and empirical solar surface models (Fontenla et al. 2006). The flow becomes more 
asymmetric as the depth of a convection zone increases (i.e., the upflows are broader 
and slower), so that there is a non-zero flux of turbulent kinetic energy, and for deep 
convection zones (^ lifp, where Hp is a pressure scale height), the turbulent energy 
flux is significant relative to enthalpy flux (~ 0.1) and oppositely directed. 

5. Future Prospects 

5.1. Comparison to Stars 

These insights are being implemented into an algorithm for stellar evolutiorlj. The idea is 
to use fully 3D phenomena, found in simulations and captured by analytic theory, by pro- 
jecting them onto a ID geometry, as used in stellar models. Unlike MLT, the algorithm is 
nonlocal and time-dependent, not static, ft should be applicable to deep, nearly adiabatic 
convection without modification. Because it have some time dependence it should be use- 
ful for models of pulsating stars, ft will replace "overshooting" and "semi-convection" 
because it uses the bulk Richardson criterion for the extent of convection (Meakin & 
Arnett (2007b)). Because the turnover flow in the convection zone is averaged over, this 
algorithm is not limited by the corresponding Courant condition, and is appropriate for 

f John Lattanzio has dubbed this the "32 ID" algorithm. 
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stellar evolution over long time scales. We emphasis that failure is possible now that free 
parameters are being eliminated, so that inadequacies of the theory will be evident. 

5.2. 3D Hydrodynamic Atmospheres 

The "321D" approach merges naturally with work on 3D atmospheres (Arnett, Meakin, 
& Young 2009 and above) and work on accretion disks (Balbus & Hawley 1998, Balbus 
2009, Blackman 2010, and Blackman & Pessah 2010). These approaches all use mean 
field equations, starting from the same general equations of mass, momentum and energy 
conservation for fluids, and use averaging to derive general properties. Because of this, 
physical processes are not introduced in patchwork fashion, but a logical necessities of the 
conservation laws. Insights into MHD in disks can spark insight into angular momentum 
transport in stars, and insights into stellar turbulence should do the same for accretion 
disk theory. As Bohdan Paczynski was fond of saying, " accretion disks are just flat stars." 

5.3. Rotation and Magnetic Fields 

Perhaps the greatest challenge for stellar evolution is the treatment of angular momentum 
transport. The rigid rotation of the Sun's radiative core, and the differential rotation of 
the convective envelope, inferred from helioseismology, seem to have been a surprise. If 
we wish to understand GRB's and hypernovae, most workers seem to assume that a key 
role is played by rotation in the gravitational collapse and explosion (an idea dating back 
to Fred Hoyle, at least). We expect to have little success if we extrapolate from the sun, 
using algorithms that give the wrong qualitative behavior. The von Neumann proposal, 
generalized to include rotation and magnetic fields, offers hope. 

Figure 1 shows the results of a first step toward understanding that problem. Our 
convection in a box simulation is continued, but with the box being rotated around the 
polar axis. The initial rotation is rigid body, so that the specific angular momentum is 
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Figure 1. Specific angular momentum versus radius. The convection zone readjusts to 
constant specific angular momentum, not rigid body rotation (Meakin & Arnett 2010). 
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quadratic in the radius. After a few turnover times, the results in Figure 1 is obtained, 
in which the specific angular momentum tends toward a constant in the convection zone, 
while remaining rigid body outside. Further, magnetic instabilities (MRI, etc.) seem to 
cause radiative regions to tend toward rigid body rotation, even if they initially have 
some other rotation law. A perusal of the literature suggests that in stellar evolution, the 
opposite is often assumed. 

6. Conclusion 

The von Neumann proposal of using computation and theory together seems to work 
well for stellar turbulence, and promises to be of value for the more complex problem 
which includes rotation and magnetic fields. Perhaps the best aspect of this approach 
is that it certainly will make new predictions of phenomena which hitherto have been 
essentially in the realm of observation only. 
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